Comparisons with PIOMAS sea ice thickness estimatesΒΆ
In this notebook, we compare our ICESat-2-derived sea ice thickness estimates with model-derived sea ice thickness from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS, Zhang and Rothrock, 2003)
import xarray as xr # For working with gridded climate data
import pandas as pd
import numpy as np
import itertools
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting
# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook
# Remove warnings to improve display
import warnings
warnings.filterwarnings('ignore')
Read in dataΒΆ
Below, weβll read in the data and restrict it to the Inner Arctic region. For the PIOMAS sea ice thickness, weβll set 0 thickness to nan to match the ICESat2 data; the ICESat-2 dataset defines anywhere with no sea ice as nan, whereas PIOMAS sets anywhere with no sea ice to 0m thickness.
# Read/download the data
book_ds = read_book_data()
# Restrict to the Inner Arctic domain
innerArctic = [1,2,3,4,5,6]
book_ds = book_ds.where(book_ds.region_mask.isin(innerArctic))
# Apply a filter to only seelct grid-cells where both ICESat-2 and PIOMAS have positive thickness (e.g. thickness greater than 1 cm)
book_ds = book_ds.where((book_ds["ice_thickness_int"]>0.01) & (book_ds["piomas_ice_thickness"]>0.01))
# Set 0 --> nan for PIOMAS thickness
#book_ds["piomas_ice_thickness"] = book_ds["piomas_ice_thickness"].where(book_ds["piomas_ice_thickness"]!=0, np.nan)
# Apply an additional filter where the gridded monthly sea ice concentration is greater than 90% - i.e. just focus on the pack ice.
book_ds = book_ds.where(book_ds["sea_ice_conc"]>0.9)
# Years over which to perform analysis
years = [2018,2019,2020]
# Select variables of interest
#book_ds = book_ds[["ice_thickness_int","ice_thickness_unc", "piomas_ice_thickness", "sea_ice_conc"]] # Grab data variables of interest
# Apply the thickness uncertainity as uper and lower values to simplify plotting uncertainity shading with Holoviews
book_ds['ice_thickness_upper']= book_ds.ice_thickness_int+book_ds.ice_thickness_unc
book_ds['ice_thickness_lower']= book_ds.ice_thickness_int-book_ds.ice_thickness_unc
Map monthly PIOMAS dataΒΆ
data_var = "piomas_ice_thickness"
#interactiveArcticMaps(book_ds[data_var], frame_width=500)
Map mean winter differencesΒΆ
Here, weβll compute winter means for PIOMAS and ICESat-2 sea ice thickness. Then, weβll use the interactiveArcticMaps function to generate maps of the winter means, along with a comparison map showing the gridcell differences.
is2_winter_means = compute_gridcell_winter_means(book_ds["ice_thickness_int"], years=years, start_month="Nov")
pio_winter_means = compute_gridcell_winter_means(book_ds["piomas_ice_thickness"], years=years, start_month="Nov")
for i in range(len(years)):
data1 = is2_winter_means.isel(time=i)
data2 = pio_winter_means.isel(time=i)
pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=data1.time.values.item())
pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title=data2.time.values.item())
diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="gridcell difference", title="Gridcell difference")
data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
display(data_pl)